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ABSTRACT 

Isothermal and self-gravitating systems bound by non-conducting and conducting 
walls are known to be unstable if the density contrast between the center and the 
boundary exceeds critical values. We investigate the equilibrium and dynamical evo- 
lution of isothermal and self-gravitating system embedded in potential well, which 
can be the situation of many astrophysical objects such as the central parts of the 
galaxies, or clusters of galaxies with potential dominated by dark matter, but is still 
limited to the case where the potential well is fixed during the evolution. As the ratio 
between the depth of surrounding potential well and potential of embedded system 
becomes large, the potential well becomes effectively the same boundary condition as 
conducting wall, which behaves like a thermal heat bath. We also use the direct A^- 
body simulation code, NB0DY6 to simulate the dynamical evolution of stellar system 
embedded in potential wells and propose the equilibrium models for this system. In 
deep potential well, which is analogous to the heat bath with high temperature, the 
embedded self-gravitating system is dynamically hot, and loosely bound or can be 
unbound since the kinetic energy increases due to the heating by the potential well. 
On the other hand, the system undergoes core collapse by self-gravity when potential 
well is shallow. Binary heating can stop the collapse and leads to the expansion, but 
the evolution is very slow because the potential as a heat bath can absorb the energy 
generated by the binaries. The system can be regarded as quasi-static. Density and 
velocity dispersion profiles from the A-body simulations in the final quasi-equilibrium 
state are similar to our equilibrium models assumed to be in thermal equilibrium with 
the potential well. 
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1 INTRODUCTION 

Thermodynamics of self-gravitating syste m is interesting 
subjec t but not clearly established yet. As IPadmanabhanI 
l|l990l ) mentioned, it is probably only fair to say that we do 
not have a systematic understanding of the self-gravitating 
system at a level similar to the kinetic theory of plasmas. Al- 
though some results from statistical mechanics may not be 
direclty applied to systems with long-range forces like grav- 
ity yBinncy & Trc maine 2008, ) , self-gravitating systems are 
correctly described by standard statistical mechanics pro- 
vided that the thermodynamica l limit is correc tly defined 
l|Padmanabhanlll990l : lKatj|2003l : IChavanidlioO^ ). 

There have been numerous attempts to understand 
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the thermodynamical behavior of self-gravitatin g system. 
Among them, related with stellar dynamics, Antonovl 
(196^) studied the entropy of self-gravitating, isothermal 
gaseous system surrounded by adiabatic rigid wall (i.e. non- 
conducting wall) and found that when the central concen- 
tration exceeds the critical value, the system can not have 
a local maximum value of entropy and leads to runaway in- 
stability. This is called A ntonov's problem sugg e sting core 



rge stmg c 
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collapse of stellar system. iLvnden-Bell fc Woo 
tended Antonov's problem to various boundary conditions 
and named the gravothermal catastrophe for this instability, 
which has been confirmed by m a ny analytical and numerica l 
works jHorowitz fc Kat3 Il978l: iHachisu fc Sugimotd [197^ : 

3 



Inagakil Il980l: ICohnI Il980l : iLvnden-Bell fc EggletonI ll98C 
Joshi et al.ll2000l ) 



The stability of isotherm al self-gravit ating system was 
first rigoro usly investigated by|Kat j l| 19781 ) , and later recon- 
sidered by IPadmanabhanI (| 19891 ) using the second variation 
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of entropy. These analyses were done in the microcanonical 
ensemble where the energy of t he system is cons erved (i.e. 
non-conducting wall). Recently IChavani i (|2002al l extended 
the work of Padmanabhan to the canonical ensemble where 
the temperature of the system is fixed (i.e. conducting wall), 
using the second variation of the free energy. 

Since the canonical distribution cannot be derived from 
the microcan onical distribution in the presence of long-range 
interactions (jPadmanabhanl Il990l ) , mean field theory has 
been used to study the thermodynamics of self-gravitating 
systems. In this perspective, self-gravitating isothermal sys- 
tem is stable only if the system is in a local maximum 
of an appropriate thermodynamical potential (i.e. the en- 
tropy in the microcanonical ensemble an d the free energy i n 
the cano nical ensemble) as m e ntioned bv lChavanid (l2002al ). 
However Ide Vega fc Sanched (l2002al lbh found the 'dilute' 
thermodynamic limit (particle number — >■ oo and volume 
V — >■ cxD, keeping N/V^^^ constant) where energy, entropy, 
the free energy are extensive. Their works justify the previ- 
ous analyses and specify the range of validity. 

Previous studies of thermodynamical description of 
self-gravitating system have mostly considered the sys- 
tems enclosed by rigid wall to prevent particle evapora- 
tion, which can be justified if a quasi-stationary condi- 
tion is satisfied (i.e. partic le evaporation rate is small). 
IVelazauez fc Guzmanl ((20031 ) replaced this rigid wall by tidal 
energy prescription and investigated Antonov's problem in 
alternative point of view, which naturally determines the 
size of the system in addition to confirming main features of 
the isothermal sphere model (i.e. core collapse and negative 
heat capacity). 

Here we propose another natural boundary condition: a 
potential well which keeps particles from evaporating. The 
completely isolated system is hard to find in astronomy and 
this type of boundary condition is often seen in different as- 
tronomical scale: for example, cluster galaxies embedded in 
dark matter potential well and dense stellar system in galac- 
tic nuclei surrounded by much larger bulge. However the sys- 
tem embedded in potential well and the effect of potential 
well to the evolution of central self-gravitating system have 
not been studied in thermodynamical point of view. There- 
fore in this work, by introducing simple model to describe 
the self-gravitating system in potential well, we attempt to 
answer the following questions: what the role of potential 
well is, how the embedded system evolves dynamically and 
what the equilibrium configuration for this embedded self- 
gravitating system can be. 

This paper is organized as follows. In Section 2, the 
previous works on the self-gravitating isothermal sphere sur- 
rounded by spherical rigid wall are briefly overviewed. We 
assign this separate section to describe the derivation of 
some formulae and provide interpretations of the previously 
known results because more detailed description of these 
previous works is necessary to describe our work which re- 
places the rigid wall by a potential well. Then in Section 3, 
we consider the potential well as a new boundary and study 
the role of potential well. In Section 4, we numerically sim- 
ulate the dynamical evolution of self-gravitating system in 
potential well. The equilibrium models are presented in Sec- 
tion 5. In the last section, the results are summarized, and 
implications and limitations of this study are discussed. 



2 SELF-GRAVITATING ISOTHERMAL 
SPHERE ENCLOSED BY 

NON-CONDUCTING AND CONDUCTING 
WALL: OVERVIEV^ 

Previous analyses on the self-gravitating isothermal sphere 
surrounded by thermally con ducting and n o n-conductin g 
wall are studied in detail bv^lPadm anabhanl (Il989l . Il990l '). 
Katj (120031). andlChayanisl ll2002a. 200^) and summarized in 



Binnev fc Tremaine I (|2008l ) . A system composed of A'^ parti- 



cles can be rep resented by one particle dis tribution function 
/ = f{x,p,t) (|Binnev fc Tremaine llioosi '). With the defini- 
tion of Boltzmann-Gibbs entropy: S = — J fXnfdPxd^p, the 
solution with extreme entropy is well known to be a spher- 
ically symmetric isothermal sphere (|Padmanabhanl 1 19891 . 
Il990f ). If we introduce the len gth, mass and energy scale 
following IPadmanabhanI (| 19891 '): 

Lo = (47rGpc/3)''", 

<^o = = 



Mo = 'inpcLo, 
kgT _ GMq 
m Lo ' 
and use new dimensionless variables: 
p(r) _ M(r) 



Lo 



(1) 



7A = /3(0-0o), (2) 



(3) 



these dimensionless variables satisfy the following relations 

, m , 2 / mn 
V = m =n^ , n = — 

Hereafter, ' symbol means the derivative with respect to ^ 
Using these relations, we obtain the Lane-Emden equation 



with the boundary condition V'(O) 
mology invariants: 

ip' m ' 



V^'(O) 



(4) 

0. Using ho- 

(5) 
(6) 



the equation describing isothermal sphere is transformed to 
the following coupled differential equations (|Padmanabhanl 
•1989). 



1 du 
u d^ 
1 dv 
vdi 



1(3 



v) 



(n-1). 



(7) 
(8) 

(9) 



^ = -5/3at(u,«) = (3,0) 
i^'iO) = 0, we obtain a spiraling 



We combine these equations to get 
u dv u ~ 1 

V du u + V ~ 3 

If we solve this equation numerically using the boundary 
conditions of n = at u 
corresponding to ^/'(O) 
curve on 11 — f plane as shown in Fig. [T] (also see Fig. 2 in 
IPadmanabhanI (|l989h 'l. The equilibrium isothermal sphere 
must exist on this curve in u ~ v plane. Since the enclosed 
mass within r of singular isothermal sphere diverges as r 
increases, it is more physically meaningful to consider the 
cut-off at radius R. Two simple boundary conditions have 
been considered: non-conducting spherical wall where no en- 
ergy is transfered and conducting spherical wall with fixed 
temperature T, where energy is transfered. 
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Figure 1. The solution of isothermal sphere in u — v plane (thick 
solid line). As the curve spirals into the point (u, v) = (1, 2), the 
density contrast between the center and the boundary increases 
and the sphere becomes singular isothermal sphere p ~ 
Two thin solid lines correspond to the energy and the temper- 
ature bound respectively. All equilibrium isothermal sphere with 
RE/GM^ = A surrounded by non-conducting wall must be lo- 
cated at the intersection between the spiral curve and a straight 
line with slope A~^. GM(5/R = rj for all equilibrium isothermal 
sphere surrounded by conducting wall has to intersect with the 
spiral curve. 



For non-conducting wall, the energy of self-gravitating 
system is conserved. Using Eqs. [5] and [6l the dimensionless 
energy A is defined as 

RE 1 / 3n 
2' 



A 



— {uo 



(10) 



where E is total energy and M is total enclosed mass within 
R. Subscript indicates the value at ^ = i^/Lo- If we rewrite 
this equation in slightly different form, we have a linear line 
on 11 — w plane with the slope given by i 



1 , 3s 
a("°"2) 



(11) 



In order for the system surrounded by non-conducting wall 
to be in isothe rmal equilibrium A must be greater than 
A. = -0.335 (lAntonovl [19621 : iBinnev fc Tremaine I l2008l : 
|Padmanabh"anlll989l . Il9m 

For conducting wall, the temperature of the system is 
conserved. The dimensionless inverse temperature r; is de- 
fined as 
GM/3 



n ■ 



R 



(12) 



Therefore, as shown in Fig. [Jl -q must be smaller than 

I ic = 2.52 for this system to be in isothermal equilibrium 
Binnev fc Tr emaine 2008; Chavanis 2002j). 

If we introduce the concept of temperature T of a stellar 
system with A'' stars, the heat capacity of the system is 

^ dE 3 - _ , 

^=dr=-2^'^- 
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Figure 2. The r)(A) curve for isothermal sphere. The heat capac- 
ity is positive until the system reaches the point Ci and becomes 
negative after Ci. After Ci, the isothermal sphere surrounded by 
a conducting wall is unstable. After C2, the isothermal sphere 
surrounded by a non-conducting wall is unstable. 



We see that the heat capacity is negative: the more energy 
the system loses, the hotter the system becomes. This ap- 
parently paradoxical phenomenon is seen not only in stellar 
system but also in any fin ite system governed by gravity 
(|Binnev fc Tremaine I [200^ ). The heat capacity C can also 
be written using A and ri jchavanis 2002a ) 

dE _ dE dp _ NkB 2 dE 
^ dT ^ 'd^dT ^ M ' lip 
Using EQS. I7l8ll0l and[T2l we get 
d\ 
di 



-NkBTj' 



:dX 
dr) 



2v^ 



Thus 

C = - 



= 7(--l)- 



-Nk, 



Au^ + 2uv - 11m -I- 3 



(14) 

(15) 
(16) 

(17) 



< 0. 



(13) 



2(u-l) 

Using A and r], we can regard the stellar system as a ther- 
modynamical system with total energy E and temperature 
T. 

We sho w ?7 as a function of A in Fig. [2] (also see 
Figure 7.1 in iBinnev fc Tremaine I l|2008l ')'). As the isother- 
mal sphere follows the solid spiral curve from the lower 
right corner, the density contrast between the center and 
the boundary D = -^^j^ increases and the heat capacity 

characterized by ^ varies over the range between —00 
to +00. If the isothermal sphere surrounded by a non- 
conducting wall passes the point C2 where D — 709, the 
system is unstable although the heat cap acity is positive 
(Antonov 1962; Bin nev fc Tremaine II2OO8I: [C havanis 2002"" 
^Horowitz & Katz 1978': 'Lv nden-Bell fc Wood ,1968, : iKat 
1978; Padmanabhan 1989, 1993). This instability originally 
introduced by lAntonovl (| 19621 ) is later called the gravother- 
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Figure 3. Upper panel: ?j as a funetion of the density contrast 
D = _f g,^ for an isothermal sphere surrounded by a conduct- 



ing wall. All points beyon d the first positi on where the = 
are known to be unstable l lChavani£ll2002al ). Lower panel: A as a 



function of density contrast D 



for an isothermal sphere 



surrounded by a non-conducting wall. All points b eyond the fi rst 
position where ^jr = are known to be unstable jKat j|l978i '). 



mal catastrophe (|Lvndeii- Bell fc Woodll 19681 ). If the isother- 
mal sphere surrounded by a conducting wall with fixed tem- 
perature T passe s the point Ci where D = 32.1, the sys- 
tem is unstable (.Binnev fc Tremaine 200j; IChavani 3 l2002al : 



iHorowitz fc Katzll978l : lKatj|l978l) However it does not leac 
to 'core-halo' structure ( Chavanig 2002a ) contra ry to the 
case of non-conducting wall ( Padmanabhaniri990l ) . 

In Fig. [3] we show the A and 77 as a function of D. Mak- 
ing Boltzmann-Gibbs entropy extreme (i.e. SS — 0) in the 
microcanonical ensemble, the equilibrium solution of self- 
gravitating system is isothermal. At every point on X{D) 
curve in the lower panel of Fi g. [3l (56' = and at critica l 
points where = , S^S = ilL\mden-Bell fc Woqdl[l9& 



|Padmanabh"arjTl989l . Il990l ). iPadmanabhanI (|l989l . Il990l l 

showed that the entropy of the system is in a local max- 
imum at any point on the branch OA (i.e. 5 S < 0) and 
the entropy of the system is in a local minimum along the 
branch AB (i.e. S^S > 0), except at A which is a saddle 
point. The branch OA is st able and the branch AB is un- 
stable (|Padmanabhaniri989l ) . It is known that al l points be- 
yond A are unstable based on previous analyses (|Kat j|l9^ . 
1 19791 ). Similarly Chavanis (2002a) showed that using the free 
energy (F — E — TS) instead of entropy S, the system in 
the canonical ensemble (conducting wall) becomes unstable 
after the first position (i.e. marked by A) where ^ = in 
the upper panel of Fig. [S] 



3 SELF-GRAVITATING ISOTHERMAL 
SPHERE EMBEDDED IN POTENTIAL 
WELL 

The boundary condition surrounding isothermal sphere 
in previous works has been either conducting or non- 
conducting wall. These boundary conditions were necessary 
to make the problem simple. Here we replace the boundary 
conditions with the realistic potential function and demon- 
strate that the potential well is similar to the heat bath if 
the potential well is deep compared with the potential depth 
of central embedded system. 

In this work, a spherical stellar system is considered to 
be embedded in potential well whose center coincides that of 
the stellar system under consideration. Then we can follow 
the same procedure of Section 2 by adding extra terms. If 
we include the external potential well (j>ext, we can write the 
potential and kinetic energy of the system as follows 



U 
K 



R 



GM{r) dM 



'0 

3M _ GM^ 3 
2^ ^ Lo 2 



dr 



dr + 



2 Z"^" dm 



dM , 
dr 

GMi 3 



(18) 



Co 



(19) 



where R is the size of embedded system and ^0 = R/Lo- 
We use Plummer potential as (p^^t and approximate it as a 
harmonic potential near the center by Taylor expansion: 



GMe 



1 + 



GAL , _ 



(20) 



where Me and a are the mass and the scale length of 
Plummer potential. Here, the harmonic approximation is 
valid as long as the embedding potential has the spatial 
scale much larger than the scale of the stellar system, as 
is often the case with the stellar syst ems in Galac t ic cen - 
ter (i.e. < .4pc core radius, see e.g. lEckart et al. I (|2005l ): 
iFieer et all (|l999l )) embedded in much larger bulge (i.e. 
« 0.56 Kpc characte ristic sca le length modeled by Dehn en 
profile iDehnenll 19931 ). see e.g. lBinnev fc Merrifieldl l| 19991 )). 
or the galaxies in cluster (i.e. radial distribution of galaxies 
modeled by broken power-law h as a scale radius 20% smalle r 
than the cluster radius, see e.g. van der Marel et al. I l|2000l ): 



iMo. van den Bosch fc White! (f201Gl )). 

Then we can rewrite Eg. 1 181 using Eqs. Ill2l and[51as 



U 



gmS 



GMoM, 



Co 



mn^dS^ 



Co 



(21) 



Then the total energy E is 



E = U + K 



GMI 



GMoMe 



GMl 
2Lo 



Co 



- d^ 



Co 



2a 



Co 
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+ 



Me Ll ^ 4 



(22) 



By defining the mean density within r &s pr = ^t'^, 



r'-' Mo 



i— = in. We can write the dimensionless 



energy A of the system as follows 
RE _ Co 



A 



GM2 



(2no^o -3mo-2Amo+Bmo^o - 



'0 



where A = ^^,B 



2Bn 
15 mjj 



)(23) 



Using u and v, this can be rewritten as 



mo 



(1.5 + A) 



1 Pe(0) 
6 p(ii) 



where pe(0) 



J_ Pe(0) _£jj_ 2 
15 p(ii) p(fl) "O 

Recalling that u 



(24) 



in, -7§T = — , we can rewrite Eg. 1241 as 



yo = 



up - (1.5 + ^) 

10 PR 



A 



(25) 



This new relation gives a linear line which is different from 
Eq. [Til In Eq. [25l the term ^ = is the ratio of 

depth of the external potential well and isothermal sphere 
at the center. Also the previous slope determined as A~^ is 

modified to (A - ^ ^|^) ~ for a given mean density of the 
embedded isothermal sphere pn and central density of the 
external potential pe(0). Since we are interested in the case 
where the potential well is deep and not much affected by 
the evolution of the central isothermal sphere, A ^ 1 and 



the slope (A 



1 Pe(0)' 



should be small (or A 



1 Pe(0) 



10 pa ) 10 PR 

should be large) but negative in order for the line to intersect 

with u — V curve. On the other hand, i^iIHi 

' PR 



_ Me h3 



' PR Mo °- 

if the size of the sphere R is assumed to be much smaller 
than the characteristic scale a of external potential and as 
a result, the term ^(-f)^ is significantly smaller than 1, 
which is appropriate assumption. Therefore A — £|M ^ \ 
and Eq. [25] approximately becomes 



^o = i{ito-(^+A)}. 



(26) 



For A S> 1, we can determine maximum A, which turns out 
to be the tangent line to the spiral curve on it — u plane in 
Fig. [2] Fig. [4] shows the same spiral curve as Fig. [2] with 
several lines determined by Eq. [26] with different values of 
A. 

For each line, there is an associated critical Ac. How- 
ever the Ac is now upper boundary in contrast to the previ- 
ous case of isothermal sphere surrounded by non-conducting 
wall, which gives lower boundary Ac = —0.335. If A is larger 
than Ac, the self-gravitating isothermal sphere surrounded 
by the potential well can not be in equilibrium. As the po- 
tential depth of surrounding potential well becomes deeper 
(larger A), the tangential line becomes close to the line 
r\c = 2.52, which sets the minimum temperature of equilib- 
rium isothermal sphere surrounded by thermally conducting 
wall. In other words, when the potential well is very deep 
compared to the potential depth of central self-gravitating 
system, the potential well behaves like a heat bath, which 
can heat up the embedded system. 
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Figure 4. The solution of isothermal sphere vau — v plane (thick 
solid line) and lines determined by Eq. 1261 for given A. Three 
thin solid lines have three different slopes determined by Ac when 
considering the external potential wells. If the value A is large, 
the ratio between the depth of external potential and potential 
of central embedded system is also large. As the external poten- 
tial depth becomes deep (large A), the line associated with the 
energy bound is close to the straight line v = rjc, which is temper- 
ature bound for isothermal sphere surrounded by heat bath. The 
inset figure magnifies the region of u — v space over the domain 
[0.9,1.2] ® [2.42,2.58] and shows that as A increases, the point 
where the tangent line and the curve meet becomes close to the 
point {u, v) = (1.0, 2.52) where the straight line v = r)c meets the 
curve. 



The heat capacity of the system surrounded by the po- 
tential well can also be written as 



C = -Nk, 



Au^ + 2uv - (11 -I- 2A)u + 3 -f 2A 
' 2(«- 1) 



(27) 



in the same way to obtain Eq. [TT] We can get the relation 
between A and r] with different A as shown in Fig. [5] Please 
note that lower left panel of Fig.[S]is the same as Fig. [5] (i.e. 
A = 0.0). There are critical points where C = or C ^ oo in 
Fig. [5] When there is no potential well {A = 0.0) or shallow 
potential well [A = 1.2), a critical point appears first at C — >■ 
cxD. If potential well becomes deep and A finally exceeds 1.5, 
a critical point first appears when C = 0. As the potential 
well becomes even deeper, the location where C = in A — 77 
plane is close to the location where C — >■ cxd. This means that 
the instability of the self-gravitating system surrounded by 
the potential well appears on the nearly same point where 
the instability of the system surrounded by heat bath occurs. 
As also seen in Fig. [4] if A is larger or smaller than 1.5, there 
is upper or lower bound value of Ac for equilibrium. 

The heat capacity of the system is shown in Fig. [5] 
for different A. Without external potential well {A — 0), 
the heat capacity becomes negative after the density con- 
trast D is greater than 32.1. Then the negative heat capac- 
ity makes the system gravitationally unstable and collapse. 
When D > 709, the heat capacity becomes positive but the 
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Figure 5. The r]{X) curve for the isothermal sphere embedded 
in the external potential well. If A < 1.5, the heat capacity is 
initially positive and grows to infinity = 0, then becomes 
negative. However if j4 > 1.5, heat capacity is initially negative 



drj 







and increases to positive value, then becomes infinity ^ 
changing its sign. And if A = 1.5, C is when the system is 
homogeneous (i.e. density contrast is 1), then grows to infinity as 
the system goes inhomogeneous density structure (i.e. increasing 
density contrast). 
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Figure 6. Heat capacity as a function of density contrast D = 
p^j^j . The heat capacities of isothermal spheres with potential 
wells diverge at the same point D = 32.1. However the points with 
C = below which the heat capacity is negative and gravothermal 
catastrophe occurs are different. As the potential well becomes 
deep (i.e. A ^ 1.5), the value of D where C = 0, is close to the 
value of D at C — > oo. 
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Figure 7. Upper left panel shows r; as a function of density con- 
trast D which is the same as the upper panel of Fig. |3] Other 
panels show A with different A, as a function of density contrast 
D. Please note that the lower panel of Fig. [3] is the case where 
A = 0. Upper right panel (A = 1.5) shows that, in contrast to the 
case of A = 0.0, the point where = 0.0 first occurs at D = 1.0 
and then appears again at D ~ 2000. Lower right panel (A = 5.0) 
shows that A reaches a upper bound beyond which the isothermal 
equilibrium state does not exist, at Z) ~ 20.0 where = 0.0. 
Lower left panel (A = 20) shows that the upper bound A appears 
at £) ~ 29.0, which is close to D = 32.1 where the instability of 
self-gravitating isothermal sphere enclosed by heat bath occurs. 
If A keeps increasing, the point where the maximum A appears 
becomes close to Z) = 32.1. 



system is known to be unstable (|Kat j [l978l . 1 19791 ). When 
the external potential becomes deep (e.g. A = 20 in Fig.|6]), 
the heat capacity of the system is negative even if it has 
homogeneous matter distribution with small D. However in 
this case, the gravothermal collapse is restrained since the 
system is embedded in deep potential well, which behaves 
like a heat bath with high temperature and heats up the 
system. Then the heat capacity becomes positive when D is 
very close to 32.1, which is the critical density contrast where 
the isothermal sphere surrounded by the heat bath becomes 
unstable. This also supports the argument that the external 
potential behaves effectively like a heat bath. 

fn Fig. [71 we show the A(D) curves of the embedded 
isothermal sphere with different A (i.e. different potential 
depth) and compare them to the 77(D) curve of the isother- 
mal sphere. The upper left panel shows the 77 as a function 
of D which is the same as the upper panel of Fig. H Other 
panels show how the A(D) of isothermal sphere embedded in 
potential well changes with different potential depth. When 
A < 1.5, there is a lower bound A for equilibrium, however, 
if A > 1.5, there is a upper bound for equilibrium. As A 
becomes large, the density contrast at upper bound A ap- 
proaches to 32.1, beyond which the self-gravitating isother- 
mal sphere within conducting wall, becomes unstable. This 
indicates that as potential well becomes deep, it approaches 
to the same boundary condition as the conducting wall (i.e. 
heat bath). 
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Table 1. Properties of initial models 



Model 


N 


Me 


i?2 




K 


Mc a 
R M 


Modell 


10000 


100.0 


10 


-28.6 


1.666 


18.62 


Model2 


10000 


100.0 


50 


-14.0 


0.450 


8.33 


Models 


10000 


100.0 


75 


-11.6 


0.369 


6.80 


Model4 


10000 


100.0 


100 


-10.1 


0.333 


5.89 



t^tot = Ws 



We 



-K 



4 DYNAMICAL EVOLUTION OF STELLAR 
SYSTEM WITHIN POTENTIAL WELL 

As shown in Section 3, the external potential well, if it 
is deep enough compared with the self-gravitating system 
in the center, increases the temperature of self-gravitating 
system and becomes effectively the same boundary condi- 
tion as conducting wall. Therefore it is interesting to study 
how self-gravitating system in potential well evolves. We use 
GPU version o f the direct A'^-body simulation code NBODY6 
l|Aarsethlll999h and simulate the dynamical evolution of stel- 
lar system embedded in the external potential well which is 
assumed to be fixed during the evolution. 

We generate four simple models. In each model, a stellar 
system following Plummer density profile is enclosed by the 
external potential wells with different depths. The central 
stellar system is generated from Plummer initial condition 
built in NB0DY6, using 10000 particles, and has the follow- 
ing form of potential-density pair. 



GM 

a 
3M 



1 + 

a 



-1/2 



-5/2 



(28) 
(29) 



The Hermite scheme in NB0DY6 require s force vectors 
to be differentiated up to the third order lAarsethl [l985l . 
l2003h . The external Plummer potential {4>e) is already im- 
plemented in NBODY6. The external Plummer potential 
has the following potential-density pair with different mass 
Me and scale length R. 



GMe 

R 

3Me 



-1/2 



-5/2 



(30) 
(31) 



471 

Please note that the notation a and R used in this section 
are difTerent from those in Section 3. 

We use the following general units used in A'^-body sim- 
ulation: 



time : 



(-4S)3/2 



length 



GM^ 

-4E 



mass : M 



(32) 



wher e G = 1,M = 1 and E = -j llHeggie fc Mathieul 
19861'). Using these units, a = 37i/16 (|Aarseth et al.lll974 



Heggie fc Hu^l2003h . Initial conditions of our models are 



listed in Table [T] Masses for external potential well are all 
set to 100 in A''-body simulation unit. Square of scale length 
R^ is different and shown in Af-body simulation unit. 

In Fig. [S] we show the evolution of Lagrangian radii 




Figure 8. The comparison between the evolution of Lagrangian 
radii of each model and isolated Plummer model (black line). 
Clearly we see that core collapse time of the system embedded in 
the potential well is longer than that of isolated Plummer model. 
When the external potential depth is deep, this discrepancy in- 
creases. 



which contain 1,5,50 and 85% of the total mass of embedded 
stellar system. The black line in each figure is the result of 
isolated Plummer model and shown for comparison. Color 
lines are the result from 4 different models. Time is scaled 
by initial half mass relaxation time (jSpitzeijl 19871 ). 



Trh.i — 



(33) 



15.4G2mplnA 

Using p ~ and N = ^, Trh i used in NB0DY6 can 

h 

be rewritten as 

Sn N <v'^ >^/^ Rl 



Trh.i — 



3 15.4G2Af2 1nA 



(34) 



where InA = ln(77V) is Coulomb logarithm determined b y 
two body relaxation and 7 is usually 0.4 (|Aarsoth 200^). 
For the isolated Plummer model, the inner Lagrangian radii 
(1,5%) decrease and the outer Lagrangian radii (50,85%) 
increase due to the gravothermal catastrophe which leads 
to the core collapse occurring at (15 ~ 16Trh,i) as seen in 
Fig. [8] Then later the inner Lagrangian radii increase due 
to the outward heat flow by the two body interaction. 

However the evolution of Lagrangian radii of the stellar 
system is retarded if it is surrounded by the external poten- 
tial well, as shown with color in each panel of Fig. [S] The 
potential depth ratio between the external potential and the 
embedded stellar system is and increases from 5.89 for 

Model4 to 18.62 for Modell, as seen in Table U NBODY6 
time step is scaled by Trh,i for each model using Eq. 1341 
As the potential well becomes deeper, core collapse time be- 
comes longer than that of isolated system or core collapse 
does not occur. 

Another interesting point is the concentration of the 
embedded stellar system compared to the isolated system. 
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Figure 9. The evolution of kinetic, potential and total energies 
of embedded stellar system for each model. Here the external po- 
tential energy is not included in the potential energy. It is easily 
seen that Modell is gravitationally unbound. Other models are 
bound by gravity and have negative heat capacity, thus expect 
to have core collapse. However the external potential increases 
the velocity dispersion of embedded system, which makes the re- 
laxation time longer than that of isolated Plummer model. Also, 
after a long time, the models suppose to reach the thermal equi- 
librium due to the interaction with the external potential well 
that behaves like a heat bath. 
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Figure 10. The snapshots of radial profile of density and velocity 
dispersion of Model2 and Model4. Model2 quickly reaches the 
thermal equilibrium with the external potential well and thus its 
velocity dispersion becomes isothermal. The final radial density 
profile of Model2 in the lower left panel is slightly changed from 
the initial density profile. Model4 experiences core collapse at t = 
SOT^h.i- Although the gradient of velocity dispersion is reduced 
with time, there is still a gradient observed in outer region (r > 
rji) of the Model4 after core collapse. The radial density increases 
by two orders of magnitude at core collapse and by several tens 
at final quasi-equilibrium state. 



Since the potential well behaves like a heat bath and in- 
creases the velocity dispersion of the embedded stellar sys- 
tem, the central region of the embedded stellar system ex- 
pands due to the increased velocity dispersion while the ra- 
dius of outer boundary is fixed because the system is con- 
fined by the external potential well. In the case of deep po- 
tential well (i.e. Modell in lower left panel of Fig. [8]), it is 
clearly seen that inner Lagrangian radii of embedded system 
are respectively larger than those of isolated system as seen 
in Fig.m 

We also show the evolution of kinetic {K), potential 
{W , without the external potential) and total energy [E) 
of embedded stellar system in Fig. [21 Ai'-body simulation 
time is scaled using Eq. [SD The initial E, K and W of 
isolated system are ~\, \ and — i respectively. For Model4, 
we see that the initial E is larger than — i, but less than 
0. Thus it is gravitationally bound and gravothermal catas- 
trophe occurs due to the negative heat capacity. The initial 
K and W are about 0.30 and —0.47 respectively. If we com- 
pare the initial potential and kinetic energies of the Model4 
to those of isolated model, we see that the kinetic energy 
significantly increases from 0.25 to 0.30. The total energy of 
the Model3 and Model2 are also negative, but more close to 
(i.e. less tightly bounded than Model4). And the diflerence 
of kinetic energy of our models from that of isolated Plum- 
mer model is larger than the difference of potential energy. 
Modell is gravitationally unbound (i.e. E is positive). Ki- 
netic energy of the Modell is much larger than the values of 
other models. As the potential depth becomes deeper, the 



kinetic energy increases due to the heating by the potential 
well. Thus velocity dispersion of embedded system increases 
and this leads to large Trh,i (see Eg. I34p . This means that 
the potential well makes the relaxation process slow. While 
the embedded stellar system heated by the potential well 
tends to expand due to increased velocity dispersion, the 
outer parts can not expand because the potential well con- 
fines the system. 

In Fig. IIOI we show the snapshots of density and veloc- 
ity dispersion profiles for two models (Model2 and Model4). 
Radius is scaled by initial half mass radius. In the left pan- 
els (Model2), black line corresponds to the initial density 
(lower panel) and velocity dispersion (upper panel) pro- 
files. The blue, green and red lines represent the profiles 
at T = 10, 20 and 50T,.h,i respectively. In the right panels 
(Model4), the black, blue, green and red lines are the pro- 
files at T = 0, 95, 100 and ISOTrh,^ respectively. Note that, 
for Model4, core collapse occurs at T = 95Trh,i- 

For Model2, we see that the velocity dispersion profile 
becomes isothermal after a few times of Trh,i- Core-collapse 
does not occur even if T = 50Trh,i, and the density pro- 
file varies slightly from the initial profile: central concen- 
tration decreases slightly and outer region is more sharply 
truncated. Model2 is less gravitationally bound system than 
Model4 (see Fig. [O]). Initially there is a gradient of velocity 
dispersion (i.e. the center is warmer than the outer). How- 
ever, as the stellar system embedded in a deep potential 
well dynamically evolves, the velocity dispersion becomes 
isothermal as a result of thermal equilibrium with the po- 
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tential well (effectively a heat bath) although the slight gra- 
dient is seen beyond Rh- From the density profile in lower 
left panel of Fig. [10] we can see that the embedded stellar 
system becomes less concentrated because the inner region 
heated by the potential tends to expand and outer region 
confined by the potential can not. 

On the other hand, we observe core collapse and gradi- 
ent of velocity dispersion profile for the Model4 in the right 
panels of Fig. 1101 As seen in upper left panel of Fig. [SI 
core collapse occurs at 95Trh,i and the density profile (blue) 
at core collapse has large concentration and shows 'core- 
halo' structure. Although Model4 is surrounded by external 
potential well, it is gravitationally bound system with neg- 
ative energy (see upper left panel of Fig. [9} . Therefore the 
self-gravity of the embedded system dominates its evolution. 
The velocity dispersion was not isothermal initially and the 
gradient of velocity dispersion still exists after core collapse, 
in contrast to the case of Model2 where the initial gradient 
of velocity dispersion is ironed out fast due to the heating 
by deep external potential well. 

From these A'^-body simulation results we see that the 
external potential well makes the relaxation process of em- 
bedded self-gravitating system slow by heating the system 
and retards core collapse, or prohibits core collapse if the po- 
tential well is deep enough. Also from the result of long term 
dynamical evolution of embedded stellar system, we expect 
a final quasi-equilibrium state of the system due to the ther- 
mal equilibrium with the external potential well acting as a 
heat bath. 

Similar phenomenon has been noticed from the study 
of dynamical evolution of two component ste llar systems 
with relatively large mass ratio.lLe3 l|l995l . l200ll ) studied the 
stellar system composed of ordinary solar mass stars and the 
black holes of 10 times higher mass. The black holes form 
compact subsystem through the dynamical friction in short 
time scale. As the central density increases the binaries form 
among black holes and eventually stops the core collapse. 
The evolution after the collapse is characterized by nearly 
static configuration since the larger stellar system composed 
of ordinary stars efficiently absorbs the heat generated by 
the binaries. Since the embedding system has much larger 
mass, the heating does not affect the surrounding system. 



5 EQUILIBRIUM MODELS FOR 

SELF-GRAVITATING SYSTEM EMBEDDED 
IN POTENTIAL WELL 

Motivated by the expectation of the quasi-equilibrium state 
in A'^-body simulation results, we consider equilibrium mod- 
els for the self-gravitating stellar system embedded in po- 
tential well. Although several equilibrium models of isolated 
system are known, among which are isothermal. King, Plum- 
mer model, a little attention is given to the equilibrium 
model of self-gravitating system embedded in potential well. 
Here we propose an equilibrium model of this system based 
on the argument discussed in Sections 3 and 4. 

A spherically symmetric stellar system with 
isotropic velocity d ispersion satisfies Jeans equation 
l|Binnev fc Tremaine \\200S, ). 



And self-gravitating system also satisfies Poisson equation 



dr 



_ _ d$ 
^ dr 



(36) 



If we consider the potential well and assume that it is deep 
enough to make the stellar system isothermal with the same 
temperature of external potential well as expected from the 
simulation results, Jeans equation may be rewritten as 



2 dp _ d4>s 
dr ^ dr 
Then, we obtain 



dr 



P = Poe 



(37) 



(38) 



where po is the central density of embedded system and, 0e 
and 03 are the external potential well and the potential well 
of self-gravitating system respectively. Thus, if we assume 
spherical symmetry, Poisson equation is 



— — (r — — ) = 47rGpoe ^ 
r'' dr dr 



(39) 



This becomes the equation of equilibrium isothermal sphere 
if (^e = 0. If we use Plummer potential for (j>e, 



R 



-1/2 



(40) 



(35) 



and solve the Poisson equation, we obtain the density profile 
of equilibrium model embedded in the Plummer potential 
well. For solving this equation, we use a normalized length 

( 

(-r(^f («) 

where we have two parameters to be set: central density 
Po and isothermal velocity dispersion of embedded self- 
gravitating system. These two parameters can usually be 
determined by observation. 

In Fig. 111! we show the equilibrium density profiles of 
isothermal sphere embedded in three different Plummer po- 
tential wells in Modell, Model2 and Model4 (see Table[l|. In 
the figure, thin solid line is isothermal sphere shown for com- 
parison. Thick green, blue and red solid lines are the equilib- 
rium density profiles for the system embedded in Plummer 
potentials in Modell, Model2 and Model4. 

As shown in Figs. [S] and [5] in Section 4, Model4 col- 
lapses due to gravothermal instability and expands later. 
Therefore the equilibrium density profile assumed to be in 
thermal equilibrium with Plummer potential in Model4 (red 
solid line in Fig.llip is concentrated, however the outer part 
is truncated due to the potential well confining the system, 
which is in contrast to the case of isothermal sphere (thin 
black solid line) . Modell whose equilibrium density profile is 
shown with green solid line in Fig. 1111 is gravitationally un- 
bound and the total energy is greater than 0.0 (see Fig. [9)1 
since the deep external potential heats the central system 
and increases its velocity dispersion (i.e. kinetic energy). The 
system is similar to the ideal gas. Thus core collapse would 
not occur. Since the potential well confines the outer part 
of the embedded system and increases the velocity disper- 
sion, the central region expands and outer boundary shrinks 
when it reaches the equilibrium. Model2 whose equilibrium 
density profile is shown with blue solid line in Fig. 1111 is the 
intermediate case and loosely bound by gravity as the total 
energy is less than, but close to 0.0 (see Fig. [9}. As shown 
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Figure 11. The equilibrium density profile of Modell (Ion- 
dashed), Model2 (dashed) and Model4 (dotted). The isothermal 
model (solid) are also shown for comparison. Model4 and Modell 
is the case of the deepest and the shallowest potential well in 
Table [T] If there is a potential well, it confines the embedded self- 
gravitating system and makes the density profile of the system 
truncated and deviated from isothermal sphere (thin solid line). 
As potential well becomes deep, the density profile becomes more 
steep in outer region. 

in Fig. [8] core collapse does not occur until the simulation 
stops at T = 78Trh,i, although it might occur after very 
long time. Since the potential well of Model2 is not as deep 
as that of Modell, the outer boundary is not declined as 
sharply as Modell. Also the central region does not expand 
as much as Modell because the increase of velocity disper- 
sion due to the heating by the external potential is not as 
significant as Modell. 

In Fig. [T2] we compare the equilibrium density pro- 
files from analytic model (dotted line) and A''-body model 
(solid line), for Model2 (blue) and Model4 (red). Final quasi- 
equilibrium density profile from A''-body simulation (black 
solid lines among the density profile snapshots of Model2 
and Model4 in Fig. [TU] were compared with equilibrium den- 
sity profiles of Model2 and Model4 in Fig. \TT\ All density 
profiles are normalized to the central density, and radius is 
scaled as follows. We estimate the po and the velocity dis- 
persion assumed to be isothermal, from density and velocity 
dispersion profiles (black lines) in Fig. 1101 and calculate the 

corresponding core radius of isothermal sphere ^ a-kGpo ^^'^ 
Model2 and Model4, where a is the radial velocity dispersion 
estimated from A'"-body simulation. Then we rescale radius 
of A^-body density profile by multiplying initial half-mass 
radius (i.e. recall that the radius in Fig. [TO]is scaled by ini- 
tial half-mass radius), and multiply the corresponding core 
radius. Now A'-body density profile radius is the same ^ as 
one in Fig. 1111 Although the discrepancy between analytic 
and A^-body models is seen at the center possibly due to the 
small number of particles in A^-body simulation, analytic 
and A-body model show the consistent result. 



Figure 12. The comparison between the equilibrium density pro- 
files from analytic model (dotted line) and A^-body model (solid 
line), for Model2 (blue) and Model4 (red). All density profiles are 
normalized to the central density, and radius is properly scaled as 
explained in the text. Although the equilibrium density profiles 
from analytic and A-body model do not perfectly agree especially 
at the center, the overall shapes are similar. 

6 DISCUSSION 

Using simple models we studied the effect of surrounding po- 
tential well to the surrounded self-gravitating system, sim- 
ulated the dynamical evolution of the system and proposed 
equilibrium models. In the following, we summarize and dis- 
cuss the result. 

6.1 The role of potential well surrounding 
isothermal self-gravitating system 

We approximate Plummer potential to harmonic potential 
near the center by Taylor expansion and investigate its ef- 
fect to the isothermal self-gravitating system at the center of 
the potential well. As the external potential becomes deeper 
compared with that of the embedded self-gravitating system, 
the external potential behaves like a conducting rigid wall 
which permits the heat exchange and conserves the tem- 
perature of self-gravitating system. If the potential depth 
ratio A is smaller than 1.5, there is a minimum dimension- 
less energy Ac in A — plane, below which the system has 
no equilibrium condition. However if A is greater than 1.5, 
there is a maximum Ac, beyond which the system can not be 
in equilibrium. As the potential depth ratio becomes large, 
the density contrast D = at A = Ac becomes close to 
32.1, which is the value for isothermal sphere enclosed by 
conducting wall with dimensionless temperature r\c = 2.52. 

Thermodynamical description of self-gravitating system 
is useful for understanding the global evolution of the sys- 
tem. Recently similar works in Se ction 2 are done using 
a general functional jTsaiii^' 198^ w hich gives the poly- 
trope with i ndex n l|Cha vanis 2002allbl; iTaruva fc Sakagamil 
l2002l , l2003al lb[). Isothermal sphere and Plummer model cor- 
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respond to the polytrope with n — >■ oo and n = 5 respec- 
tively. However, the maximization of Tsallis functional at 
fixed mass and energy is a condition of dynamical stability 
rather than thermodynamical stability (,Chavan is 2004a). In 
this context, polytropic distribution is justified as a particu- 
lar steady solution of the coUisionless Boltzmann equation. 
Furthermore, in this dynamical interpretation, Tsallis func- 
tional is not an entropy. 

Strictly speaking, the self-gravitating system does not 
have the thermodynamic limit where usually particle num- 
ber — » oo and vo lume V ^ oo, keeping N/V constant 
de Vega fc Sanchez 2002a'). However, de Vega & Sanchez 



de Vega fc Sanchej|2002al ibl) found the 'dilute' thermody- 
namic limit (particle number N ^ oo and volume V — >■ cx3, 
keeping N/V^^^ constant) where energy, entropy, the free 
energy are extensive. This study provides a justification 
of taking thermodynamical approach to describe the self- 
gravitating system, which is useful to understand important 
physics using much less expensive computational resource 
than numerical simulation. 



6.2 The effect of potential well to the dynamical 
evolution of the embedded self-gravitating 
system 

We generate self-gravitating stellar system and surround it 
using external Plummer potential. NBODY6 simulates the 
evolution of the system and shows the consistent results with 
the argument in Section 3. The potential well retards the 
relaxation process by heating the embedded stellar system 
and increasing its velocity dispersion. Thus if the embedded 
system has a positive energy E, it behaves like an ideal gas 
and does not experience the gravothermal catastrophe. On 
the other hand, the system with negative energy eventually 
experiences core collapse by gravothermal catastrophe, al- 
though core collapse time of the system is larger than that 
of isolated stellar system. It is because, as the kinetic en- 
ergy of stars interacting with the potential increases, the 
system becomes loosely gravitationally bound and Trh,i of 
the embedded system increases (see Eg. I34p . 

The evolution of Lagrangian radii of our model shows 
that the deep external potential makes the embedded system 
gravitationally unbound and core collapse does not occur 
as seen in Fig. [S] From the energy exchange between the 
embedded system and surrounding potential well as seen 
in Fig. [§1 it is suggested that the potential well heats the 
embedded system and increases its kinetic energy. Therefore 
the total energy of embedded system can be positive if the 
potential well is deep enough. 

From the time evolution of density and velocity disper- 
sion profiles of the embedded system in the external poten- 
tial well, it is very likely that the embedded system is in 
thermal equilibrium with the potential well. In deep poten- 
tial well, we show that the velocity dispersion profile be- 
comes isothermal due to the heating by the potential well. 
The inner part of density profile becomes less concentrate 
and the outer part becomes steeper than the isolated sys- 
tem. In the shallow potential well, we see that core collapse 
occurs, and the velocity dispersion profile is not isothermal 
and has a gradient. 

These simulations are based on simple model where the 
possible interaction between central system and surrounding 



potential well is not considered. In order to understand how 
the embedded system co-evolves with potential well, we need 
to implement the potential well using large number of stellar 
particle instead of fixed potential function. However, for the 
deep potential well, our simple approach with fixed potential 
would be a good approximation. 



6.3 Equilibrium configuration for self-gravitating 
system embedded in potential well 

Based on the conclusion that if the potential well is a heat 
bath, the embedded self-gravitating system is in isothermal 
equilibrium with the potential well, we derived the equilib- 
rium density profiles of self-gravitating system embedded in 
potential well by solving Jeans equation and Poisson equa- 
tion. These equilibrium density profiles are similar to A'^- 
body simulation results. 

Especially these equilibrium models are physically mo- 
tivated by radial distribution profile of galaxies in c luster, 
which is often modelled by King profile (|King|[l966l ). King 
model is a good fit for distribution of galaxies near the clus- 
ter core. However King model looses stars by tidal energy 
cut-off. This picture is unrealistic in the case of galaxies in 
cluster, which are embedded in the deep potential well of 
cluster dark matter halo. Galaxies in cluster are bound and 
hard to escape from cluster potential well. Our model can be 
more realistic description for galaxy clusters. While the King 
profile drops at the outer part due to tidal energy cut-off, 
our model profile drops because the potential well confines 
the embedded system and keeps stars from escaping out. 
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